Fermi surface of a trapped dipolar Fermi gas 
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Under the framework of the semi-classical theory, we investigate the equilibrium state properties 
of a spin polarized dipolar Fermi gas through full numerical calculation. We show that the Fermi 
surfaces in both real and momentum spaces are stretched along the attractive direction of dipolar 
interaction. We further verify that the deformed Fermi surfaces can be well approximated by 
ellipsoids. In addition, the deformation parameters slightly depend on the local real and momentum 
space densities. We also study the interaction strength dependence of the energy and real and 
momentum space densities. By comparing them with variational results, we find that the ellipsoidal 
ansatz usually generates accurate results for weak dipolar interaction; while under strong dipolar 
interaction limit, notable discrepancy can be observed. Finally, we map out the stability boundary 
of the system. 
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I. INTRODUCTION 

The novel long-range and anisotropic character of 
dipole-dipole interaction has stimulated significant in- 
terest on studying dipolar Bose- Einstein condensate [l|. 
Experimentally, the first dipolar condensate was realized 
in Cr atoms, which possess a large magnetic dipole mo- 
ment ■ Utilizing Feshbach resonance to tune the scat- 
tering length to near zero, the dipolar effect was also ob- 
served in condensates of Rb Q and Li 0] atoms. For spin 
polarized fermionic atoms, even though the contact inter- 
action vanishes, the Pauli exclusion makes the magnetic 
dipole-dipole interaction difficult to observe. However, 
the experimental progress on trapping and cooling po- 
lar molecules, in particular the recent success in making 
high phase-space-density fermionic KRb gas 0, 0] , opens 
up a new avenue for realizing strongly correlated states. 
Owing to the anisotropic dipolar interaction, we expect 
more novel features to be added to the quantum Fermi 
gases. 

On the theoretical work of dipolar Fermi gases, You 
and Marinescu [7[ first pointed out that the p-wave paired 
BCS states could be achieved for fermionic atoms inside 
an external dc field. Subsequently, the critical temper- 
ature of the superfluid transition and its relation to the 
trap geometry were investigated [1, . Other theoretical 
work also includes studying the strongly correlated states 
in rapidly rotating trap [nj, HH and the possible biaxial 
nematic phases [12j |. 

For dipolar Fermi gas in normal phase, the equilibrium 
state properties [lU , low-lying collective excitations [3| , 
and expansion dynamics [lj| have been studied based on 
the semi-classical theory. In those works, it was assumed 
that the momentum distribution is isotropic, which im- 
plicitly eliminates the possibility to explore the effect of 
exchange dipole-dipole interaction. Miyakawa et al. (l6j 
then proposed an ellipsoidal variational ansatz for the 
phase space distribution function, which contains a pa- 
rameter characterizing the deformation of the Fermi sur- 



face in momentum space. It was shown that the exchange 
dipolar interaction induced the deformation in momen- 
tum distribution and destabilized the dipolar Fermi gas. 
Following this work, the collective excitation and free ex- 
pansion of a dipolar Fermi gas was also studied (l7j . 

In the present work, we investigate the ground state 
properties of a trapped dipolar Fermi gas in normal 
phase. We justify the use of the ellipsoidal variational 
ansatz by numerically finding the phase space distribu- 
tion function. We show that, in both real and momentum 
spaces, the Fermi surfaces can be ap pro ximated by ellip- 
soid as proposed by Miyakawa et al. flq ]. In addition, we 
find that the deformation parameters are weakly depen- 
dent on the spatial and momentum coordinates. 

This paper is organized as follows. In Sec. [Til we in- 
troduce our model and briefly outline the semi-classical 
theory for ultra-cold Fermi gas. In Sec. IIII1 we present 
the numerical algorithm employed in this work. In Sec. 
IVi we find numerically the phase space distribution func- 
tion of a trapped dipolar Fermi gas, from which we study 
the equilibrium state properties of the system and com- 
pare them with those obtained variationally. Finally, we 
conclude in Sec. IVl 



II. FORMULATION 

We consider a system of N spin polarized fermionic 
particles with permanent dipole moment d at zero tem- 
perature. For simplicity, we assume that all dipoles are 
polarized along z-axis by an external field. The dipole- 
dipole interaction potential becomes 
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where 9 T . 
vector r 



./ is the angle between positive z-axis and the 
r' and Cd — rjd 2 / (Atteo) with rj being a pa- 
rameter continuously tunable within the range [— \ , 1] . 
The parameter r\ can be realized using a fast rotating 
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orienting field [l8| . It not only makes the sign of dipolar 
interaction changeable, it can also be used to completely 
switch off the dipolar interaction if ip equals to 54.7°, 
the 'magic angle'. For spin polarized system, the contact 
interaction corresponding to s-wave scattering vanishes. 

In second quantized form, the Hamiltonian of the sys- 
tem takes the form 



H 



h 2 V 2 
2m 



U ho (r) 



$(r) 



+- J drdr'ft(f)ft(r')V d (T - r')^(r')^(r), (2) 

where m is the mass of the molecule, ip(v) is the field 
operator, and U^ (r) = \m(w\x 2 + <jj\y 2 + uj 2 z 2 ) is the 
external trapping potential which is assumed to be axi- 
ally symmetric. In our numerical simulation, to make the 
volume of trapping potential to be a constant, we shall fix 
the geometric average of trap frequencies Hi = (ui^uiz) 1 ' 3 . 
To characterize the shape of the trapping potential, we 
define the trap aspect ratio as A = w 2 /wi, such that the 
harmonic trap now becomes 



LWr) = ^ 2 A- 2 /V + AV), 



(3) 



where p 2 — x 2 + y 2 . 

Under Hartree-Fock mean field theory [19], the wave 
function of the system is a Slater determinant |$) = 
c i c 2 ' ' ' c jv |0) w hh c\ being the creation operator gen- 
crating a molecule at orbital <j>i(r). Using single-particle 
reduced density matrix p(r,r') = ($| ip^ ( r ')V'( r ) the 
total energy E = ($\H\$) can be expressed as 

E = j dr|^-[V r -V r ^(r,r')] r=r , + p(r, r)f/ ho (r)| 
+~ J drdr'p(r,r)p(r',r')V d (r-r') 
~ J drdv'p(r',v)p(T,r')V d (r-r'). 
To proceed further, we shall employ the Wigner function 



/(r,k) = J dse- ik »p(r+~,r-^), 



(4) 



which is the Fourier transform of the single-particle re- 
duced density matrix. Under classical limit, the total 
energy can now be expressed as as functional of Wigner 
function 



E[f] = 



1 
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drdk 



h 2 k 2 
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l — J drdkdr'dk7(r,k)/(r',k')U d (r-r') 



2(2tt) 



drdkdk7(r,k)/(r,k')Vd(k - k') 



(5) 



where Vd(k) = Cd^-{i cos 2 6^ — 1) is the Fourier trans- 
form of the dipolar interaction with 9^ being the polar 
angle of k. The two terms in first line of Eq. ([5]) corre- 
spond to, respectively, the kinetic and potential energies, 
the second line represents the direct dipolar interaction 
energy, and the last line is the exchange dipolar interac- 
tion energy. 

The Wigner function can be interpreted as the phase 
space distribution function, which is normalized to the 
total number of molecules 



N = 



1 



(2tt) 3 



drdk/(r,k). 



(6) 



Integrating over k and r, we obtain, respectively, the real 
space density n(r) = [l/(27r) 3 ]J dk/(r, k) and the mo- 
mentum space density fi(k) = [l/(27r) 3 ] J rfr/(r, k). At 
zero temperature, the phase space distribution function 
of Fermi gas takes the form of Heaviside step function 



/(r,k) = e(Mr,k)), 



(7) 



where Sf(y, k) = defines a closed surface which we refer 
to as the Fermi surface in phase space. For those phase 
points enclosed by the Fermi surface, we have /(r, k) = 1, 
otherwise, /(r, k) = 0. 

For the simplest case, we may straightforwardly adopt 
the local density approximation which assumes that the 
Fermi energy only depends on the local real space density 
n(r), i.e., 



4 sph) (r,k)^ 2 n(r)] 



,2/3 



(8) 



Due to the spherical symmetry of the momentum space 
distribution, the exchange dipolar vanishes [131 ] . such 
that the total energy becomes a functional of n(r) only 



Jj(sph) r 
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+ £7 ho (r)n(r) 



-- J drdr' n(r)n(r')Vd{r — r'). 



(9) 



The real space density can now be obtained either varia- 
tionally or numerically by minimizing 

An improvement over Eq. ([8]) was made by Miyakawa 
et al. [la ] , who proposed an ellipsoidal ansatz in momen- 
tum space 



4 ell) (r,k)^ 2 n(r)] 
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k„ + a kt 
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— , , k 2 and the dimensionless parameter 



where k 2 — k 2 

a characterizes the deformation of Fermi surface in mo- 
mentum space. Clearly, the Fermi surface in momentum 
space forms an ellipsoid, 



3l 



(11) 
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where K± = a 1 / 2 (67r 2 n) 1/3 and K\\ = a' 1 (6n 2 n) 1 ^ 3 . 
Moreover, the real space density is assumed to have an 
inverted parabolic shape 



( 48 ^) 1/3 ir-©vv+/3v) 



3/2 



(12) 



where a — yj h/(muj) is the length of the harmonic os- 
cillator. The dimensionless variational parameters (3 and 
7 correspond to, respectively, the deformation and com- 
pression in real space. The phase space distribution is 
then completely characterized by a, f3, and 7 which can 
be determined by minimizing the total energy. The im- 
mediate consequence of allowing the deformation in mo- 
mentum space is that the exchange dipolar interaction 
energy is nonzero, which has a direct impact on the sta- 
bility of the system 

Furthermore, based on the variational ansatz Eqs. (jlOp 
and (|12|. the Fermi surface in real space takes the form 



P 

R 2 



z 



1 



(13) 



where R± = /3 1 / 2 [3n/(4 7 r)] 1 / 3 andi?,, = /T 1 [3n/(4 7 r)] 1 / 3 . 
It is worth pointing out that, in the noninteracting limit, 
the real space density takes the exact form 



n 



1 

6^ 



(48JV) 



1/3 



(A- 2 V + A 4/3 ; 



-1 3/2 



(14) 



Comparing it with Eq. P^|) . we see that, in the absence 
of interaction, the real space deformation /3q — A 2 / 3 is 
completely determined by the trap aspect ratio. 

In the present paper, we shall determine the phase 
space distribution function of a trapped dipolar Fermi 
gas using full numerical method and compare our results 
with the variational results using ellipsoidal ansatz. Be- 
fore doing that, let us first briefly outline the numerical 
algorithm employed in this work. 



III. NUMERICAL ALGORITHM 

Due to the cylindrical symmetry of the system, the 6- 
dimensional phase space distribution function reduces to 
a 4-dimensional one as 



/(r,k) = f(p,z,k p ,k z ). 



(15) 



To obtain /, we shall numerically minimize the free en- 
ergy 



T[f\ = E[f] - fiN, 



(16) 



where the chemical potential /1 is introduced to fix the 
total number of particles. The minimization is carried 
out using simulated annealing method. 



To proceed, we introduce the 4-dimensional grids of 
extent [0,R] x [-Z,Z] x [0,K p ] x [-K Z ,K Z ] with total 
number of grid points N r x N z x N r x N z . The values 
of R, Z, K p , and K z are chosen such that we know for 
sure f{p,z,k p ,k z ) — for phase space point outside the 
grid extent. Since in both real and momentum space the 
system is symmetric with respect to xy-plane, the grid 
extent is reduced to [0, R] x [0, Z] x [0, K p ] x [0, K z ], For 
simplicity, we usually set N r = N z in our numerical cal- 
culation, and typically, they are within the range of 40 to 
80. For the simulated annealing method, we also intro- 
duce a fictitious temperature T which controls the speed 
of convergence. Now, the simulation can be implemented 
straightforwardly as follows: 

1. For an initial phase-space distribution /, a start- 
ing temperature T, and a chemical potential /i, we 
calculate the free energy F[f] ■ 

2. To generate a new phase space distribution func- 
tion /', we select a target phase space point gi = 
(pi,Zi,k p i,k z i) and make a trial move by setting 
f( gi ) = 1 (0) if /(«?,) = (1). The free energy 
difference /S.J- = T[f] — J-[f] is then calculated. If 
AJ 7 < 0, the trial move is accepted. Otherwise, we 
accept it with probability e~~ A:F / T . 

3. Step 2 is repeated under the fixed temperature T 
until the number of trial moves reaches an upper 
bound -ZVtri, or the number of accepted moves 
reaches another upper bound Nacc- Generally, 
Nacc is much smaller than -/Vtri- In case either 
of these two conditions is satisfied, we lower the 
temperature T by a small fraction, say 10%. 

4. Step 2 and 3 are repeated until the free energy con- 
verges. 

To expediate the convergence of the simulation, we al- 
ways choose the target grid points for trial moves at the 
vicinity of the Fermi surface in step 2, as for a grid point 
gi deep inside the Fermi surface, it is unlikely that f{gt) 
will change during the simulation. Accordingly, the typ- 
ical values of Nacc an d -ZVtri are proportional to the 
number of grid points on Fermi surface. Finally, we re- 
mark that the chemical potential fi remains unchanged 
once it is selected at the beginning of our numerical sim- 
ulation. As a consequence, the total number of particles 
N can only be calculated after we obtain the final phase 
space distribution function. 



IV. RESULTS 

To present our results, it is convenient to rescale all 
quantities into dimensionless forms. To this end, we in- 
troduce the dimensionless units as follows: N 1 / 6 ?! for 
length, 27r7V 1 / 6 a~ 1 for wavevector, and N 4 ^ 3 hu) for en- 
ergy. Now, the phase space distribution function is nor- 
malized to unit, J drdkf(r, k) = 1. The real and momen- 
tum space densities are expressed as n(r) = J dkf(r, k) 
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A 



sionless quantities. The dimenionless kinetic and po- 
tential energies become 2tt 2 J drdkk 2 /(r, k) and E pot — 
hj dr\~ 2 / 3 (p 2 + X 2 z 2 )n(r), which can be calculated 
straightforwardly using numerical integration. Defining 
a dimcnsionless dipolar interaction strength 

a 3 huj ' 

the direct dipole-dipole interaction energy can be ex- 
pressed as 



FIG. 1: The trap aspect ratio dependence of the critical 
dipolar interaction strengths from numerical calculation (solid 
lines). As a comparison, we also plot the variational results 
(dashed lines) [la ]. 

and n(k) = [l/(2n) 3 ]J dr/(r,k), respectively. Here and 
henceforth, we adopt the same notation for the dimen- 



E dir = — J dvdr n(r)n(r ) — — 3 , 

which involves only the real space density. In cylindrical 
coordinate, E^jr can be evaluated efficiently using Han- 
kel transform [2(|. Finally, the dimensionless exchange 
dipolar interaction energy takes the form 



2 " d J drdkdk7(r,k)/(r,k')(l -3cos 2 6> k _ k ,) 

dpdzdkpdk z dk' p dk' z pk p k' p f{p, z, k p , k z )f(p, z, k' p , k' z ) 



3 

(2^) 4 £ d 



3 f k% k 



!Y2 



^(k 2 - k' 2 ) 2 + 2(k 2 + k' 2 )(k z - k' z ) 2 + (k z - k' z y 



where to obtain last two lines, we have analytically inte- 
grated out the azimuthal variables of r, k, and k'. At first 
sight, this 6-dimensional integral may look formidable for 
numerical integration, as one can only tackle it using bru- 
tal force. Using the fact that, in each step of our numeri- 
cal simulation, we only update the value of f(p, z, k p , k z ) 
on a single phase space point, the numerical integration 
can be done efficiently. 



After formulating our problem into the dimensionless 
form, it becomes clear that the control parameters in our 
system reduce to the trap aspect ratio A and dipolar in- 
teraction strength Ed- While the dependence of all phys- 
ical quantities on N can be easily found by converting 
them back to the dimensional forms [l3j] . We investigate 
below the A and Ed dependences of the equilibrium prop- 
erties of our system. In particular, we shall compare our 
results with those obtained variationally and justify the 
validity of the variational approach. 



A. Stability 

Due to the partially attractive character of dipo- 
lar interaction, there exist critical dipolar interaction 
strengths beyond which the system becomes unstable. 
This fact has been pointed out in several previous stud- 
ies [l3|, EE EH ■ In Fig. [TJ we present the stability dia- 
gram of a dipolar Fermi gas on the X-£d parameter plane 
based on full numerical calculations. Compared to the 
variational results [161 ] . the critical dipolar interaction 
strength is significantly lower. The reason that varia- 
tional approach fails to predict the stability boundary 
is because the simple ellipsoidal ansatz, Eqs. (|10p and 
(fT2|) . is incapable to capture the local collapse induced 
by strong dipolar interaction. 

To illustrate the structure of local collapse, we plot 
in Fig. [5] the typical intermediate real and momentum 
space densities from our numerical simulation. The cor- 
responding parameters are A = 10 and Ed = 2.6 which fall 
into the unstable region. Unlike n(r) for a stable config- 
uration [Fig. [3] (a) and (c)] which is a smooth function of 
the spatial coordinate, the real space density shown here 
oscillates violently. Corresponding to the sharp peaks in 
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FIG. 2: The intermediate results for the real (a) and momen- 
tum (b) space densities of an unstable system with A = 10 
and Ed = 2.6. The total energy diverges eventually as we 
continuously carry out the simulation. 

n(r), the momentum space density possesses a very long 
small valued tail, as compared to n(k) of a stable system 
[Fig- El (b) and (d)]. We emphasize that the result shown 
in Fig. [5] does not represent converged real and momen- 
tum space densities, the corresponding total energy of 
the system diverges eventually as we continuously carry 
out the simulation. 

Finally, as we have properly taken into account the 
exchange interaction energy, which trends to destabilize 
the system, the critical dipolar interaction strengths pre- 
sented in Fig. [T] are also lower than those predicated nu- 
merically by using the spherical ansatz Eq. ([8]) [H, EEl • 

B. Real and momentum space densities 

The real space density can be observed directly when 
the in-situ detection is available for the system. We 
present, in Fig. [3] (a) and (c), the typical behaviors of 
real space densities corresponding to different control pa- 
rameters A and Ed- For a cigar-shaped (pancake-shaped) 
trap, the peak real space density increases (decreases) as 
one increases the dipolar interaction strength. This in- 
dicates that the overall dipolar interaction is attractive 
(repulsive) in a cigar-shaped (pancake-shaped) trap. In 
contrast to the real space density, the peak momentum 
space density [Fig. [21(b) and (d)], decreases (increases) as 
one increases the dipolar interaction strength for a cigar- 
shaped (pancake-shaped) trap. As a comparison, wc also 
plot the real and momentum space densities obtained 
through variational method. Clearly, when dipolar in- 
teraction is weak, the agreement between two methods 
is very well; while under strong dipolar interaction, the 
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FIG. 3: (color online). Real and momentum space densities 
for A = 0.1 [(a) and (b)] and 10 [(c) and (d)]. The solid and 
dashed lines correspond to, respectively, the numerical and 
variational results. For A = 0.1, in ascendent (descendent) 
order of the peak real (momentum) space density, Ed = 0.1, 
0.5, and 1. For A = 10, in descendent (ascendent) order 
of the peak real (momentum) space density, the interaction 
strengths are Ed = 0.1, 1, and 2. 
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FIG. 4: (color online). \/(kl)/{kl) as a function of dipolar 
interaction strength for various trap geometries. 

discrepancy in peak density can be as high as 10%. 

If we switch off the dipolar interaction [l8| to let the 
gas expand ballistically, the momentum space density 
can be observed directly from the timc-of-fiight image. 
Therefore, the quantity n = \/ (k%)/{k%) represents the 
deformation of the expanded cloud. We present the dipo- 
lar interaction strength dependence of n in Fig. 0] corre- 
sponding to different trap geometries. In noninteract- 
ing case or when the interaction is isotropic, the ex- 
panded cloud always has a spherical shape. However, the 
anisotropic feature of dipolar interaction always stretches 
the momentum space density along attractive direction 
of the dipolar interaction. In addition, the deformation 
of momentum space density only weakly depends on the 
trap aspect ratio. 
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FIG. 5: (color online). The direct (a) and exchange (b) dipo- 
lar interaction energies for X = 0.1 (solid lines), 1 (dashed 
lines), and 10 (dash-dotted lines). 



C. Energy 

In Fig. 03(a) and (b), we present the Sd dependence of 
direct and exchange interaction energies for various trap 
geometries. Similar to dipolar Bose-Einstein condensate, 
the direct interaction energy strongly depends on the ge- 
ometry of the trapping potential: when Ed > 0, E^ is 
positive (negative) for oblate (prolate) trap, indicating 
that the overall direct interaction is repulsive (attrac- 
tive); while in a spherical trap, -Edir is always attractive 
as a result of the stretch of the gas along the attractive 
direction of dipolar interaction. On the other hand, the 
exchange dipolar interaction is always negative, which re- 
flects the fact that the momentum distribution is always 
stretched along the attractive direction of dipolar interac- 
tion. Moreover, except for in a spherical trap, where the 
magnitude of E c ^ c is comparable to the that of -Edin the 
exchange interaction energy is usually below 30% of the 
direct interaction energy in an highly anisotropic trap. 

Now, we turn to study the release energy E re \ which is 
the sum of the kinetic and interaction energies and can 
be measured experimentally. In Fig. [6l we plot the re- 
lease energy as a function of dipolar interaction strength 
for various trap geometries. As a result of the small 
exchange interaction energy, the behavior E re \ is quite 
similar to that obtained numerically with the spherical 
ansatz ([5]) [l5| ■ The dashed lines in Fig. [H] correspond 
to the variational results. It can be seen that even the 
discrepancy increases at the strong interaction limit, it is 
still below 1%. Moreover, for the total energy, the agree- 
ment between numerical and variational methods can be 
even better. Therefore, as long as we are in the stable 
region, the variational ansatz can be safely adopted to 
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FIG. 6: (color online). The dipolar interaction strength de- 
pendence of the release energy for various trap aspect ratios. 
The solid and dashed lines correspond to, respectively, the 
numerical and variational results. 
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FIG. 7: (color online). The Fermi surfaces in momentum 
space 5V(0, 0, k p , k z ) = (a) and in real space Sf(p, z, 0, 0) = 
(b). The corresponding parameters (X,£d) are indicated in 
the figures. The dashed lines in (a) and (b) are plotted using, 
respectively, Eqs. (|17[1 and (|19p with K'± n and iijii obtained 
from the numerical results. 



calculate the energy of the system. 



D. Fermi surface 



The numerical studies of the real and momentum space 
densities and energies of a trapped dipolar Fermi gas al- 
lows us to justify the use of variational method from the 
point of view of the global physical quantities. To gain 
more insight into the validity of variational method, we 
shall construct the Fermi surface of the system in phase 
space which enables us to compare in detail the behavior 
of the variational deformation parameters (a and 0) with 
numerical results. 
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FIG. 8: (color online). The real density dependence of K'±, 
K'u, and a for A = 0.1 and ea = 0.5. The corresponding 
variational a is plotted as dotted line. 
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FIG. 9: (color online). The momentum space density depen- 
dence of R' ± , Rf„, and /3'\- 2/z for A = 10 and e d = 0.5. The 
corresponding variational /3A _2//3 is plotted as dotted line. 



1. Fermi surface in momentum space 

In Fig. [7] (a) , we plot the Fermi surfaces in momentum 
space, Sf(0, 0, k p , k z ) = 0, corresponding to various con- 
trol parameters (A, Ed)- Apparently, the Fermi surfaces in 
momentum space are anisotropic, resembling an ellipsoid 
described by equation 



K P 

K? 



Kf 



1, 



(17) 



where K' ± and K'» correspond to, respectively, K± and 
Kn in the variational ansatz Eq. (fTTj) . By reading out the 
values of K' ± and K'u from our numerical results, we plot 
Eq. (JT7J) using dashed lines in Fig. [7] (a). The remark- 
able agreement between 5(0, 0, k p , k z ) = and Eq. (fT7| 
suggests that the Fermi surface in momentum space can 
be well approximated by an ellipsoid. We emphasize that 
the control parameters used in Fig. [7] fall into the strong 
interaction region which are close the stability boundary 
(see Fig. [T}. For weaker dipolar interaction strength ed, 
the agreement can be even better. 

We also find that K' ± and K'„ only depend on spatial 
coordinate through the real space density, namely K' ± ^ 
are only functions of n(r). Furthermore, as shown in Fig. 
HI both K' ± and K'u are roughly proportional to n 1 / 3 , in 
analog to the behavior of K± and Kn. To reveal more 
details of the deformation of Fermi surface, we define the 
ratio of K± and Kii as 



K' 

v /3/2 = f\L 



K 



(18) 



Apparently, a' corresponds to the variational parame- 
ter a in Eq. (|10p which characterizes the deformation of 
Fermi surface in momentum space. The typical result 
of a' is plotted in Fig. [5] as dash-dotted line. We im- 
mediately see that the momentum space Fermi surface is 
stretched along z-axis (a' < 1), in agreement with predic- 
tion of variational method. However, in contrast to the 
assumption that a is a constant (dotted line in Fig. [5]) in 



the ellipsoidal ansatz, we find that a' is a decreasing func- 
tion of n, which indicates that higher real space density 
associates with larger momentum space deformation. 



2. Fermi surface in real space 

As shown in Fig. [7] (b) , the Fermi surface in real space 
can be similarly expressed as an ellipsoid of the form 



Rf 



Ft' 2 

U \\ 



1. 



(19) 



In addition, the shape of Sp(p, z, 0, 0) = very sensi- 
tively depends on the trap geometries. To characterize 
the deformation of Fermi surface in real space, we define 







3/2 



R'u ' 



(20) 



^From Eq. (|14[) . we immediately realize that the quantity 
0/^—2/3 re p re sents the real space Fermi surface deforma- 
tion purely induced by dipolar interaction. In Fig. [9l we 
plot the momentum space density dependence of R' ± and 
R'u . Similar to the momentum space case, we also find 

that R' ± || are roughly proportional to [n(k)] 1 / 3 . How- 
ever, the deformation parameter /?' clearly shows that 
higher momentum space density corresponds to larger de- 
formation of Fermi surface in real space. In addition, the 
dipolar interaction alway stretches the real space Fermi 
surface along z-axis (/3'A~ 2//3 < 1). We point out that 
the irregular behavior of /?' at low h end is caused by the 
limited resolution of the grids used in numerical simula- 
tion, which becomes more important when the momen- 
tum space density is small. 



V. CONCLUSIONS 

To conclude, we have studied the equilibrium state 
properties of a spin polarized dipolar Fermi gas based on 
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the semi-classical theory. Employing the simulated an- 
nealing method, we obtain numerically the phase space 
distribution function by minimizing the total energy of 
the system. We confirm that the Fermi surfaces in both 
real and momentum space can be well approximated by 
ellipsoids, which are stretched along the attractive direc- 
tion of dipolar interaction. However, in contrast to the 
ellipsoidal variational ansatz in which the deformation 
parameters are assumed to be constants, we find that 
they weakly depend on the local real and momentum 
space densities. We also study the dipolar interaction 
strength dependence of the real and momentum space 
densities. We find that the results from variational cal- 
culation agree with numerical results when the dipolar 
interaction is weak; while for strong dipolar interaction, 



notable discrepancy is observed. Finally, we map out 
the stability boundary based on numerical calculation. 
The numerical critical dipolar interaction strengths are 
significantly lower than those predicted variationally. 

We thank Liang He for helpful discussion. This work is 
supported by NSFC (Grant No. 10674141), National 973 
program (Grant No. 2006CB921205), and the "Barren" 
program of Chinese Academy of Sciences. 

Note added. During the preparation of this 
manuscript, we become aware of several work on studying 
the ground state, sound propagation, and expansion dy- 
namics of dipolar Fer mi g ases [2ll [22l l23l l2~i| . In particu- 
lar, Ronen and Bohn [21| obtain the exact Fermi surface 
of a homogeneous system through numerical calculation. 
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